Assignment pre-series04 - Python¶

By Chloé Laignel et Amélie Goubé¶

1. Generate a known signal¶

  1. write down the equation of a sinusoid at 1Hz.
  2. create a 1D time vector with 1000 samples and a sampling frequency of 100 Hz.
  3. using time, create a signal that is a 2D numpy array with the dimensions (x, y, z) over columns : use a sinusoid at 1Hz, with Pi/4 shift from x to y and to z.
  4. reproduce the output below (i.e., the text and the 2D and 3D plots)
In [1]:
conda install -c plotly plotly

2 channel Terms of Service accepted

Channels:
 - plotly
 - defaults
Platform: osx-64
Collecting package metadata (repodata.json): 
- 
\ 
| 
/ 
- 
\ 
| 
done
Solving environment: - 
\ 
| 
/ 
done

# All requested packages already installed.


Note: you may need to restart the kernel to use updated packages.
In [2]:
import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go

# -----------------------------------------
# Question 1 – Generate a known signal
# -----------------------------------------

# Sampling parameters
fs = 100              # sampling frequency (Hz)
f0 = 1.0              # signal frequency (Hz)
A = 1.0               # amplitude

# Phase shift between axes
phase_shift = 0.25 * np.pi   # 0.25π rad = 45 degrees

# Create a time vector from 0 to 3 seconds with exactly 100 samples
t = np.linspace(0, 3, 100).reshape(-1, 1)

# Build a 3D signal:
# x(t): phase = 0
# y(t): phase shifted by +phase_shift
# z(t): phase shifted by +2*phase_shift
s_x = A * np.sin(2 * np.pi * f0 * t + 0)
s_y = A * np.sin(2 * np.pi * f0 * t + phase_shift)
s_z = A * np.sin(2 * np.pi * f0 * t + 2 * phase_shift)

# Stack the components into one array (N, 3)
S = np.hstack([s_x, s_y, s_z])

# Print shapes (as requested in the assignment)
print("t.shape =", t.shape)
print("S.shape =", S.shape)
print("phase_shift =", phase_shift)

# -----------------------------------------
# 2D Plot (x, y, z over time)
# -----------------------------------------
plt.figure(figsize=(10, 5))
plt.plot(t, S[:,0], label='x')
plt.plot(t, S[:,1], label='y')
plt.plot(t, S[:,2], label='z')
plt.title("3D Sinusoidal Signal (2D projection)")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.legend()
plt.grid()
plt.show()

# -----------------------------------------
# 3D Plot (x, y, z trajectory)
# -----------------------------------------
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(S[:,0], S[:,1], S[:,2])
ax.set_title("3D Sinusoidal Trajectory")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

# limits of the axes
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-1, 1)

ax.set_xticks(np.arange(-1, 1.5, 0.5))
ax.set_yticks(np.arange(-1, 1.5, 0.5))
ax.set_zticks(np.arange(-1, 1.5, 0.5))

plt.show()

t.shape = (100, 1)
S.shape = (100, 3)
phase_shift = 0.7853981633974483
No description has been provided for this image
No description has been provided for this image

2. Estimate the derivative of the signal¶

  1. write down the math to get the derivative of a sine wave.

Given a sine function: $$ s(t) = A \sin(2 \pi f t + \phi) $$

Its analytical derivative is:

$$ \frac{ds}{dt} = 2 \pi f A \cos(2 \pi f t + \phi) $$

  1. write down the equation of the derivative of signal along each axis.

For a 3D signal:

$$ \mathbf{S}(t) = \begin{bmatrix} x(t) \\ y(t) \\ z(t) \end{bmatrix} $$

Each component is a sinusoid with a phase shift:

$$ \dot{x}(t) = 2 \pi f A \cos(2 \pi f t + \phi_x) $$

$$ \dot{y}(t) = 2 \pi f A \cos(2 \pi f t + \phi_y) $$

$$ \dot{z}(t) = 2 \pi f A \cos(2 \pi f t + \phi_z) $$

  1. create a function that estimates the time derivative of a signal:

We can estimate the derivative numerically using two methods:

1. Forward difference:

$$ \dot{s}(t_i) \approx \frac{s(t_{i+1}) - s(t_i)}{t_{i+1} - t_i} $$

2. Central difference:

$$ \dot{s}(t_i) \approx \frac{s(t_{i+1}) - s(t_{i-1})}{t_{i+1} - t_{i-1}} $$

Note: At the boundaries, forward/backward difference is used to keep the same array length.

In [3]:
# ---------------------------------------------------------
# Derivative estimation functions (Forward and Central)
# ---------------------------------------------------------

def derivative_forward(signal, time):
    """
    Estimate the derivative using the forward difference method.
    Maintains the same length as the input signal by using
    a backward difference for the last sample.
    """
    s = np.asarray(signal)
    t = np.asarray(time).reshape(-1)
    
    if s.ndim == 1:
        s = s.reshape(-1, 1)  # convert to 2D
    
    N, D = s.shape
    ds = np.zeros((N, D))
    
    # Forward difference for all points except the last
    dt = np.diff(t)
    ds[:-1, :] = (s[1:, :] - s[:-1, :]) / dt[:, None]
    
    # Backward difference at the final point
    ds[-1, :] = (s[-1, :] - s[-2, :]) / (t[-1] - t[-2])
    
    return ds


def derivative_central(signal, time):
    """
    Estimate the derivative using the central difference method.
    Uses forward/backward difference at the boundaries.
    """
    s = np.asarray(signal)
    t = np.asarray(time).reshape(-1)
    
    if s.ndim == 1:
        s = s.reshape(-1, 1)
    
    N, D = s.shape
    ds = np.zeros((N, D))
    
    # Central difference for interior points
    for i in range(1, N - 1):
        dt = t[i + 1] - t[i - 1]
        ds[i, :] = (s[i + 1, :] - s[i - 1, :]) / dt
    
    # Forward and backward differences at the boundaries
    ds[0, :] = (s[1, :] - s[0, :]) / (t[1] - t[0])
    ds[-1, :] = (s[-1, :] - s[-2, :]) / (t[-1] - t[-2])
    
    return ds
  1. apply the function to the signal and plot the result.
In [4]:
# Compute derivatives if not already done
dS_central = derivative_central(S, t)  # central difference

# Plot original signal + derivative on the same figure
plt.figure(figsize=(10,6))

# Plot original signal (x, y, z)
plt.plot(t, S[:,0], label='x (signal)', color='blue')
plt.plot(t, S[:,1], label='y (signal)', color='green')
plt.plot(t, S[:,2], label='z (signal)', color='red')

# Plot derivative of the signal (x, y, z)
plt.plot(t, dS_central[:,0], '--', label='dx/dt', color='cyan')
plt.plot(t, dS_central[:,1], '--', label='dy/dt', color='lime')
plt.plot(t, dS_central[:,2], '--', label='dz/dt', color='magenta')

# Labels and title
plt.xlabel("Time (s)")
plt.ylabel("Amplitude / Derivative")
plt.title("Signal and its derivative (superimposed)")
plt.legend()
plt.grid()
plt.show()
No description has been provided for this image

3. Estimate the tangential velocity of the signal¶

3.1 Definition

The tangential velocity is the magnitude of the instantaneous rate of change of the position vector.
In other words, it is the Euclidean norm of the derivative of the signal vector:

$$ v_\text{tan}(t) = \|\dot{\mathbf{S}}(t)\|_2 = \sqrt{\dot{x}^2(t) + \dot{y}^2(t) + \dot{z}^2(t)} $$

3.2 Formula using linear algebra

Let the signal be a matrix:

$$ S \in \mathbb{R}^{N \times d} $$

with N samples and d dimensions.
The derivative is:

$$ \dot{S} \in \mathbb{R}^{N \times d} $$

Then the tangential velocity vector is:

$$ v = \text{norm}(\dot{S}, \text{axis}=1) \in \mathbb{R}^N $$

3.3 Python function to compute tangential velocity

In [5]:
def tangential_velocity(signal, time, derivative_method='central'):
    """
    Compute the tangential velocity of a signal and return it as a 2D numpy array.
    
    Parameters:
        signal: ndarray of shape (N, d)
        time: ndarray of shape (N,) or (N,1)
        derivative_method: 'central' or 'forward'
        
    Returns:
        v: ndarray of shape (N, 1) containing the tangential velocity
    """
    s = np.asarray(signal)
    
    # Compute derivative
    if derivative_method == 'central':
        ds = derivative_central(s, time)
    elif derivative_method == 'forward':
        ds = derivative_forward(s, time)
    else:
        raise ValueError("derivative_method must be 'central' or 'forward'")
    
    # Tangential velocity = L2 norm of derivative along each row
    v = np.linalg.norm(ds, axis=1)
    
    # Return as 2D column array
    return v.reshape(-1, 1)

3.4 Test the function

In [6]:
# 1D signal
s1 = (np.sin(2*np.pi*1*t)).reshape(-1)
v1 = tangential_velocity(s1, t)
print("1D test: s1.shape =", s1.shape, "v1.shape =", v1.shape)

# 2D signal
s2 = np.hstack([s_x, s_y])
v2 = tangential_velocity(s2, t)
print("2D test: s2.shape =", s2.shape, "v2.shape =", v2.shape)

# 3D signal (already S)
v3 = tangential_velocity(S, t)
print("3D test: S.shape =", S.shape, "v3.shape =", v3.shape)

# Plot tangential velocity

# Compute tangential velocity if not already done
v3 = tangential_velocity(S, t)

# Create figure with 2 subplots (stacked vertically)
fig, axs = plt.subplots(2, 1, figsize=(10, 8), sharex=True)

# --------------------------
# Top subplot: original signal
# --------------------------
axs[0].plot(t, S[:,0], label='x', color='blue')
axs[0].plot(t, S[:,1], label='y', color='green')
axs[0].plot(t, S[:,2], label='z', color='red')
axs[0].set_ylabel("Amplitude")
axs[0].set_title("3D Signal Components")
axs[0].legend()
axs[0].grid()

# --------------------------
# Bottom subplot: tangential velocity
# --------------------------
axs[1].plot(t, v3, label='Tangential velocity', color='purple')
axs[1].set_xlabel("Time (s)")
axs[1].set_ylabel("Velocity")
axs[1].set_title("Tangential Velocity")
axs[1].legend()
axs[1].grid()

# --------------------------
# Adjust layout and show
# --------------------------
plt.tight_layout()
plt.show()

# Offset the tangential velocity to display it below the signal
offset = -1.5                  # vertical offset (adjust if needed)
v3_offset = v3.flatten() + offset  # flatten in case v3 is (N,1)

plt.figure(figsize=(10,6))


# --------------------------
# Plot the original signal (x, y, z)
# --------------------------
plt.plot(t, S[:,0], label='x (signal)', color='blue')
plt.plot(t, S[:,1], label='y (signal)', color='green')
plt.plot(t, S[:,2], label='z (signal)', color='red')

# --------------------------
# Plot the tangential velocity (offset)
# --------------------------
plt.plot(t, v3_offset, '--', label='Tangential velocity (offset)', color='purple')

# Add a horizontal reference line to indicate zero level of velocity
plt.axhline(offset, color='purple', linestyle=':', alpha=0.5)

# Labels, title, grid and legend
plt.xlabel("Time (s)")
plt.ylabel("Amplitude / Velocity")
plt.title("Signal and Tangential Velocity (offset to avoid overlap)")
plt.legend()
plt.grid()
plt.show()

1D test: s1.shape = (100,) v1.shape = (100, 1)
2D test: s2.shape = (100, 2) v2.shape = (100, 1)
3D test: S.shape = (100, 3) v3.shape = (100, 1)
No description has been provided for this image
No description has been provided for this image

4. Create an interactive plot with the plotly library¶

In [7]:
import plotly.graph_objects as go
import numpy as np

# Assume S, v3, t are already defined

# Flatten v3 in case it is (N,1)
v3_flat = v3.flatten()

# Create figure
fig = go.Figure()

# Add traces for the signal
fig.add_trace(go.Scatter(x=t.flatten(), y=S[:,0], mode='lines', name='x (signal)'))
fig.add_trace(go.Scatter(x=t.flatten(), y=S[:,1], mode='lines', name='y (signal)'))
fig.add_trace(go.Scatter(x=t.flatten(), y=S[:,2], mode='lines', name='z (signal)'))

# Add trace for tangential velocity (dashed line)
fig.add_trace(go.Scatter(x=t.flatten(), y=v3_flat, mode='lines', name='Tangential velocity',
                         line=dict(dash='dash', color='purple')))

# Update layout
fig.update_layout(
    title="Signal and Tangential Velocity (Interactive)",
    xaxis_title="Time (s)",
    yaxis_title="Amplitude / Velocity",
    legend_title="Legend",
    hovermode="x unified"
)

# Show interactive plot
fig.show()

# Export to HTML
fig.write_html("signal_and_velocity.html")